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Abstract. In this work we construct multigrid preconditioners to accelerate the solution process 
of a linear-quadratic optimal control problem constrained by the Stokes system. The first order 
optimality conditions of the control problem form a linear system (the KKT system) connecting the 
state, adjoint, and control variables. Our approach is to eliminate the state and adjoint variables 
by essentially solving two Stokes systems, and to construct efficient multigrid preconditioners for 
the Schur-complement of the block associated with the state and adjoint variables. These multigrid 
preconditioners are shown to be of optimal order with respect to the convergence properties of the 
discrete methods used to solve the Stokes system. In particular, the number of conjugate gradient 
iterations is shown to decrease as the resolution increases, a feature shared by similar multigrid 
preconditioners for elliptic constrained optimal control problems. 
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1. Introduction. Over the last decade, the computational community has shown 
a growing interest in devising fast solution methods for large-scale distributed optimal 
control problems constrained by partial differential equations (PDEs). Optimal con- 
trol problems constrained by the Stokes system form a stepping stone in the natural 
progression from the - now classical - Poisson-constrained test problem to problems 
constrained by more specialized and complex PDE systems modeling fluid flow such 
as Navier-Stokes, non-Newtonian Stokes, or the shallow water equations. Optimal 
control problems constrained by such PDE models play important roles in real-world 
applications, such as modeling of ice sheets or data assimilation for ocean flows and 
weather models. 

In this article we consider an optimal control problem with a cost functional of 
tracking- type: 

min J(u,p,f) = ^\\u- ud\\l 2( ny + lr\\P ~ Pd|[l 2 (Q) + d|/lli 2 (n)a > 
subject to the constraints 

'-Au + Vp = f in 0, 

divu = infi, (1.2) 
u = on dfl, 

where Q is a bounded jsolygonal domain in K 2 . The purpose of the control problem 
is to identify a force / that gives rise to a velocity u and/or pressure p to match a 
known target velocity Ud, respectively pressure pd- Since this problem is ill-posed, 
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we consider a standard Tikhonov regularization for the force with the regularization 
parameter /3 being a fixed positive number. The constants j u , "f p are nonnegative, 
not both zero. 

The main goal of this article is to construct and analyze optimal order multigrid 
preconditioners to be used in the solution process of (ll.ll) - (|1.2p . Over the last few 
years a significant amount of work has been devoted to developing multigrid methods 
for optimal control problems. An overview of this research and further references can 
be found in a review article by Borzi and Schulz 1 . However, fewer works are dedi- 
cated specifically to optimal control problems constrained by the Stokes system. For 
example, recently, Rees and Wathen have proposed two preconditioners for the 
optimality system in the distributed control of the Stokes system, a block-diagonal 
preconditioner for MINRES and a block-lower triangular preconditioner for a non- 
standard conjugate gradient method. We note that there are several papers in the 
literature on finite element error analysis for the optimal control of the Stokes equa- 
tions (see, e.g., jTOJ [T^l S] and the references therein). Our paper focuses on the 
solution of the linear system that arises in the discretization process, which is not 
addressed in these papers. However, for completeness, we also prove an a priori error 
estimate for the optimal control since the cost functional in (ll.lj) includes a pressure 
term, using standard techniques similar to those used in [101 1121 [5] . 

Since the cost functional in is quadratic, the KKT system is a linear saddle- 
point problem connecting the state, adjoint, and control variables. Solution methods 
for these problems typically fall into two categories: the all-at-once approach takes 
advantage of the sparsity of the system, but has the disadvantage that the matrix 
is indefinite. On the other hand, Schur-complement strategies may lead to smaller 
systems that may also be positive definite, but the sparsity is lost. Our approach is to 
eliminate the state and adjoint variables by essentially solving two Stokes systems us- 
ing specific methods (see [7]), and to construct efficient multigrid preconditioners for 
the Schur-complement of the block associated to the state and adjoint variables. The 
constructed multigrid preconditioners arc related to the ones developed by Draganescu 
and Dupont [5] , and are shown to be of optimal order with respect to the convergence 
properties of the discrete finite element methods used to solve the Stokes system. In 
particular, the number of conjugate gradient iterations is shown to decrease as the 
resolution increases, a feature shared by similar multigrid preconditioners for elliptic- 
constrained optimal control problems. One word on the optimality of the precondi- 
tioner: the usual notion of optimality, especially in the context of multigrid, refers to 
the cost of the solution process being proportional to the number of variables. We 
argue that for this problem, an unpreconditioned application of conjugate gradient 
(CG) in conjunction with an optimal multigrid solve for the Stokes system already 
satisfies this notion of optimality. In the current context, multigrid preconditioners 
actually can perform better than that, and optimality refers to the order of approx- 
imation of the operator under scrutiny - in this case the reduced Hessian - by the 
multigrid preconditioner, as shown in Theorem 13. II 

The paper is organized as follows: In Section [2] we introduce the discrete opti- 
mal control problem and prove finite element estimates that will be needed for the 
multigrid analysis. Section [3] contains the main result of the paper, Theorem 13.11 
which refers to the analysis of the two-grid preconditioner; furthermore, the extension 
to multigrid preconditioners is briefly discussed. In Section 3] we present numerical 
experiments that illustrate our theoretical results, and we formulate some conclusions 
in Section [5] 
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2. Discretization and convergence results. The strategy we adopt is to first 
discretize the optimal control problem then optimize the discrete problem. To define 
a discrete problem based on a finite element approximation of the Stokes system we 
briefly recall the usual weak formulation of the Stokes equations. Define the spaces 

v = Hl{n)\ 

M = L 2)0 (Q) =L L 2 (n) : jf p = oj, 
and the bilinear forms o:FxF->l and b : V x M — > R as 

&(<£,p) = - / pdiv<^. 
Jn 

Throughout this paper, we write (•, •) for the inner product in L 2 (f2) or L 2 (f2) 2 ac- 
cording to context, and similarly for the norms, if there is no risk of misunderstanding. 
The weak solution (u, p) G V x M of (jl.2[) is the solution of 

a(u, tp) + b(tp, p) = (/, 0) Vtp e V, 

b(u, i>) = W> e M. 

For / G i? _1 (fi) 2 the problem has a unique solution [5]. Moreover, if is a convex 
polygon and / G £ 2 (^) 2 , then u G H 2 (Q) 2 , p G ff 1 (f2) [9], and there exists C = 
C(O) > such that 

\\u\\H^ + \\V P \\ L2{n) <C\\f\\ L2m 2. (2.1) 

Throughout this paper we will assume O to be convex, so that the H 2 - regularity of 
the Stokes problem is ensured. Furthermore, the target velocity field Ud is assumed 
to be from L 2 (f2) 2 and the target pressure pd from M. 

We introduce the solution mappings U and V of the continuous state equation 
defined such that for any / G L 2 (f2) 2 the following holds: 

a(Uf,<p*) + b(0,Vf) = (f,0) and b{Uf,^)=Q V(^)eFxM. 

The mapping U, considered as a linear operator in L 2 (r2) 2 , is compact and self-adjoint, 
as 

{Ufi , /I) = a{Uh ,Uf 2 ) = (/i , uf 2 ) V/i ,/ 2 ei 2 (r») 2 . 
We denote by V* : L 2 , (^) — > L 2 (f2) 2 the adjoint operator of V, defined by 

With this notation, the problem (|l.l[) - (jl.2[> is written in reduced, unconstrained form 
as 

min J(f) = ^\\Uf- u d \\l 2m , + ^\\Vf-p d \\l 2m + §11/11,(0). ■ (2-2) 
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The Hessian operator associated to the reduced cost functional J is given by 

Hp= lu U*U + lp r*P + /3I . 

Note that the solution of the minimization problem is obtained as the solution 
of the normal equation 

Hpf = lu Wu d + lp V* Pd . (2.3) 

The goal of this paper is to design an efficient multigrid algorithm for solving the 
discrete version of (|2.3I) . 

2.1. Finite element approximation. We consider a shape regular quasi-uniform 
quadrilateral mesh 3~h of fi, and we assume that the mesh results from a coarser 
regular mesh 2?ih from one uniform refinement. We use the Taylor-Hood Q2 — Qi 
finite elements to discretize the state equation. The velocity field u is approximated 
in the space V h ° = V% R H^(rt) 2 , where 

V h = {v h £ C(Cl) 2 : vh\t G Q 2 (T) 2 for T £ ,%}, 

and the pressure p is approximated in the space 

M h = {q h e C(O) n i 2 , (O) : q h | T e Q X (T) for T £ 

where Qk(T) is the space of polynomials of degree less than or equal to k in each vari- 
able |5| . The control variable / is approximated by continuous piecewise biquadratic 
polynomial vectors in the space Vh- 

Remark 2.1. For convenience, we choose here the quadrilateral Q2 — Qi Taylor- 
Hood elements, however, our analysis can be extended to triangular P2 — Pi elements 
as well as other stable mixed finite elements. 

For a given control / E L2(0) 2 , the solution (uh,Ph) 6 K° x Mh of the discrete 
state equation is given by 

a{u hl h ) + b(ip h ,Ph) = (f,if h ) V<Ph G V°, 

b(u h ,tph) = Vtp h €M h . 

Let W/j and "P^ be the solution mappings of the discretized state equation and U'h, 
Vt their adjoints, defined analogously to the continuous counterparts. Furthermore, 
denote by nh : Li(£t) 2 — > Vh the ^-orthogonal projection onto Vh- The discretized, 
reduced optimal control problem reads 

min Jh(fh) = ^-Whfh u'}\\ 2 + ^\\V h f h -p h d f + f llAf, (2.4) 
fh£V h III 

where u^, p\ are the I/2-projections of the data onto Vh, respectively Mh- 

Let us investigate the structure of the algebraic system associated to the dis- 
cretized optimal control problem. Let {0j}*j =1 and {V'fcjfcLi be the basis functions of 
the spaces Vh and Mh, respectively. Furthermore, assume that {<£j}™ = i are the basis 
functions of V® for some n < p. We expand the discrete solutions Uh and ph as 



n m 

Uh(x) = ^2ujipj(x), p h (x) = ^p fc Vfc(a:), 
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and we approximate the control by fh{ x ) — 2j=i ^j ( Pj( x )- Let A G R™ xn and B G 
R mxn be the matrices related to the bilinear forms, A = [a^] = a(<pi, tpj), B = [6y] = 
b(0i,ipj). Furthermore consider the mass-matrix M^= [Jq </?i • $j] G R pxp and its 
submatrices M.,- = M.-(l : n, 1 : p) G K™ xp and M s = M^l : n, 1 : n) G W ixn . The 
discrete state equation takes the form 



A B T 
B 



«/ 




After discretizing, the reduced cost functional becomes 

Jh = y(u - u rf ) T M s (u - u d ) + ^(p - Pd ) T M p (p - p d ) 



^f T M f -f, 

2 / ' 



where M p = [ J„ V'iV'j] an d u^, p,j are the coefficient vectors in the expansions of 
and pf in V°, respectively M^. 
Let us introduce the matrices 



A B T 
B 



M 









7 P M P 



'ML,-' 




and Xd 



After dropping the constant terms the problem becomes 



min -f T (L T S _1 MS _1 L + /3M.M - f T L T S -1 Mx d , 
f 2 ^ 



which reduces to solving the linear system 



(I/S^MS^L + /3M f -)f = L T S _1 Mx d . 



(2.5) 



Note that the system matrix is dense, thus (I2.5[) has to be solved using iterative meth- 
ods, and for increased efficiency we need high-quality preconditioners. In Section O 
we construct and analyze a multigrid prcconditioner for the system 



(/3I + M^ 1 L T S- 1 MS- 1 L)f = M--, 1 L T S- 1 Mx d , 
which is obtained from (|2.5p by left-multiplying with ML 1 . 



(2.6) 



We should remark that the system (|2.6[) can be obtained from the KKT system 
associated with the discrete constrained optimal control problem associated to 
(|1.2[) by block-eliminating the velocity, pressure, and the Lagrange multipliers. So (12.61) 
is in fact a reduced KKT system. 

2.2. Estimates and convergence results. We first list some results on the 
Stokes equations and their numerical approximation which are needed for the multi- 
grid analysis. 

Theorem 2.2. There exist constants C\ = C\{U, SI) and C2 = CafP, SI) such 
that the following hold: 
(a) smoothing: 



l|W/1l<Ci||/W-W V/eL 2 (!]) 2 , 
IIP/1 <C 2 ||/V W V/eL 2 (!l) 2 ; 



(2.7) 
(2.8) 



() 
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(b) approximation: 

\\Uf-U h f\\ <C x h 2 \\f\\ V/eI 2 (J]f, (2.9) 

and 

\\Vf-r h f\\<C 2 h\\f\\ V/6i 2 (flf; (2.10) 

(c) stability: 

\\U h f\\ <Ci || f\\ and \\V h f\\< C 2 || f\\ V/eL 2 (!l) 2 . (2.11) 

Proof. 

(a) The inequality (|2.7p is a straightforward consequence of (|2.1[) (see [5j Corol- 
lary 6.2]), while (|2 .8[) follows immediately from Brezzi's splitting theorem, 
see, e.g., HJ Theorem 4.3]. 

(b) This is a standard approximation result, see, e.g., |13j . 

(c) Follows immediately from the estimates in (a) and (b). 

□ 

We state without proof the following well-known result. 

Lemma 2.3. The following approximation of the identity by the projection holds: 

- n h )f\\ H - k(n) , < Ch k \\f\\ V/ £ L 2 (ft) 2 , (2.12) 

for k = 0, 1, 2, 3, with C constant independent of h. 

Remark 2.4. Theorem \2.2\ and Lemma \2.3\ imvlv that there are constants C\ — 
Ci(fi,W) and C 2 = C 2 (Q,V) such that 

\\U{I--K h )f\\<C x h 2 \\f\\ V/eI 2 (f]) 2 (2.13) 

and 

\\V(I-n h )f\\<C 2 h\\f\\ V/eL 2 (U) 2 . (2.14) 

Lemma 2.5. Let tt^ : L 2 q(Q) — > Mh be the L 2 -orthogonal projection onto Mh- 
Then 

\((I-n h )q,Vf)\<Ch\\q\\\\f\\ Vg G L 2 , (O),/ e I 2 (0) 2 , (2.15) 

with C constant independent of h. Proof. We have 

|((J - n h )q, Vf)\ = |((7 - Z h )q,Vf- V h f)\ <]](/- TT h )q\\ \\Vf- V h f\\ 
OTTot 

< C%]]]]/|| VgeL 2 , (!l),/el 2 (J]) 2 . □ 

Denote 

G u = U*U , G p = V*V , G h u = W h U h , Gp = V* h V h . 



Lemma 2.6. The following approximation properties hold: 
h h {G h u-Gu)f\\ <Ch 2 11/11 VfeV h , 



(2.16) 
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and 

hh(Gp ~ G P )f\\ < Ch\\f\\ VfeV h , (2.17) 

for some constant C independent of h. 
Proof. For / G Vh we have 

\(n h (Gi-G u )f,f)\ = \\\U h ff-]\Uff\ 

< \\U h f-Uf\\(\\U h f\\ + \\Uf\\) 
<Ch 2 \\f\\\ 

and (|2.16l) follows from the symmetry of 7T/j(G^ — G u ). 
Similarly, we have 

\(n h (G>; - G p )fJ)\ = \{G h p fJ) v - (G p fJ) v \ 

= \(V h f,V h f) M -(Vf,Vf) M \ 

= \\\nf\\ 2 -\\vf\\ 2 \ 

< \\V h f-Vf\\(\\V h f\\ + \\Pf\\) 
<Ch\\f\\ 2 Vfev h , 

and (|2.17l) follows from the symmetry of Hh{G p — G p ). □ 
From the definition of Hp follows 

/3||/1| 2 <(^/,/)<(/3 + C)||/|| 2 V/eI 2 (fi) 2 , (2.18) 

with C = C(U,V,fl), and a similar estimate holds for the discrete Hessian 

ff£ = /3/ + 7 u G£+ 7p G£, (2.19) 

which shows that the condition number of Ha is bounded uniformly with respect to 
h, but potentially increasing with (3 J. 0. The inequality (|2.18[) also implies that the 
cost functional J is strictly convex and has a unique minimizer given by 

/""" = H^{ lu Wu d + lp V* Pd ). 

Similarly, the minimizer of the discrete quadratic is 

fr n = {B h v)-\iuUiu h d + lv p* hV h d ). 

In the following theorem, we show that /™ n approximates f mm to optimal order in 
the L2-norm. 

Theorem 2.7. There exists a constant C = C(fl,U,V) independent of h such 
that for h < h ((3, fi, U, V) we have the following stability and error estimates: 

< ll/ mln H + j(luh 2 \\u d \\+ lp h\\p d \\), (2.20) 

u/r n - r in \\ < ^(7^ 2 (ii^n + \\f mm \\)+ipH\\ Pd \\ + \\f mm \\)). (2.2i) 
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Proof. Let e h = f™ n - f mm . We have 

HpS h = 031 + 7 . U G. U + j P G p )fr i ~ H f mm 

= 08J + lu G u + lp G p )fr n ~ luU*u d - lp V* Pd 

= 08/ + 7u G£ + lp G h p )fr n + 7«(G U - G^)/r" 1 + 7p(G p - G h p )fr n 

- -f u U*u d - l P V*p d 

= 7< ,Z^tZ5 + 7p 7>^5 + lu {G u - G h u )fr i + J P (G P - G h p )fr n 

- -f u U*u d - l P V*p d 

= 7u [(G„ - G h u )fr n + (K - U*)u h d -W{I- 7r h )u d ] 

" v ' 

+ lp [(G p - G h p )flr + (K - V*)p d -V*(I- n h )p d ] , 



where we used that u d = ■nuUd and p h d = irtPd- Furthermore, 



and 



P\\e h \\ < (Hpe h , e h ) < J U (A U , e h ) + j p (A p , e h ), 



f2~9ll . lf2TT3l 

(A u: e h ) < Ch 2 \\u d \\\\e h \\ + ((G U -G*)fr n ,?h) 

d-^)e h ±v h ^ii^im^ii + {7Th{Gu _ G l)f^^ h e h ) 

+ (GJr n ,(I-7T h )e h ) 

< Ch 2 \\e h \\(\\u d \\ + ||/r n ||) + {Ufr n M(I - 7r h )e h ) 

^ Ch 2 \\e h \\(\\u d \\ + \\fr n \\), 

with C = C(U, il), where we have also used the fact that U is self-adjoint. Similarly, 
we get 

(A p , e h ) = (p h d , (V h - V)e h ) - ((/ - n h ) Pd , Ve h ) + ((G p - G h p )fr n ,e h ) 



< Ch\\p d \\\\e h \\ + ((G P -G h p )fr n ,e h ) 
(l-, h )e h ±v h CHpdmh] \ + {7Th{Gp _ Gfrjfr**,^) 
+ (G p fr n ,(I-n h )e h ) 

™ ch\\e h \\(\\ Pd \\ + \\fr n \\) + {vrr\ni - 



j2TT4t 

< ch\\s h \\{\\ Pd \\ + \\fr n \\\ 

with C = C(V, ri). These estimates yield 

C 



i/r n - r m \\ < -p(iutf{\\u d \\ + n/r»n) + 7 ^omi + ii/ri). 
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Since 



u/ri < n/ mm n + ii/r"-/ 

c 



< \\r m \\ + j(iuh 2 (\\u d \\ + \\fr n \\) + imwpa + \\ft in \\), 

we obtain ([2720]) and (ET21) for h sufficiently small. □ 

3. Two-grid and multigrid preconditioner for the discrete Hessian. In 

this section, we use the multigrid techniques developed in [3] to construct and analyze 
a two-level symmetric preconditioner for the discrete Hessian Hp defined in (|2.19[) . 
The extension to a multigrid preconditioner follows the same strategy as in [5] and is 
further explained in great detail in [BJ. 

For the remainder of this paper we consider on Vh the Hilbert-space structure 
inherited from L 2 (^l) 2 ■ Furthermore, we consider the Z/2-orthogonal decomposition 
Vh = V~2h®W2h and let n 2 h be the L 2 -projector onto V 2 h- The analysis in [5] suggests 
that Hp is well approximated by 

T| <^ H} h n 2h +f3(I-7r 2h ). 

For completeness we briefly recall the heuristics leading to the definition of Tt. As 

usual in the multigrid literature, for fh £ Vh we regard ir 2 hfh as the "smooth" com- 
ponent of fh, and (/ — iT2h)fh as the "rough" or "oscillatory" component; so the 
projector (/ — iz 2 h) extracts the "oscillatory" part of a function in Vh- If we write 
Hp = Hq + (31 and take into account the "smoothing" properties of Hq (these are 
due to the compactness of the operator Hq which Hq approximates), it follows that 
the products (/ — iT 2 h)HQ and Hq (I — TT 2 h) are almost negligible. So 



H% = (7r 2fc + (/ - 7r 2h ))(H$ + (3I)(ir 2h + -K2h)) 

w 7r 2h (flf + PIfah + 13(1 - K 2h ) . (3.1) 

Furthermore, when applied to the "smooth" component n 2 hfh 01 a function //, , it is 
expected that HQir 2 hfh ~ H ir 2 hfh ~ H^x^fh, hence the idea to replace in (|3.ip 
Hq by Hq H , which gives rise to Tp. 

Since n 2 h is a projection, (T^) _1 is computed explicitly as 

d ^ f {T^ = {Hf)-\ 2h +^{I-^ h ). 

We propose Lp G £(Vh) as a two- level preconditioner for Hp. To assess the quality 
of the preconditioner we use the spectral distance introduced in [5], defined for two 
symmetric positive definite operators T\,T 2 e C(Vh) as 



d h (T u T 2 ) = sup 

wev h \{o} 



(Txw,w) 



(T 2 w, w) 



(3.2) 



If Lp is a preconditioner for Hp then the spectral radius p(I — LpHp), which is an 
accepted quality-measure for a preconditioner, is controlled by the spectral distance 
between Lp and (Hp)^ 1 (see Lemma IA.2I in Appendix lAl for a precise formulation). 
The advantage of using the spectral distance over p(I — LpHp) is that the former is 
a true distance function. 
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Theorem 3.1. For h < h ((3,Cl,U,V) there exists a constant C = C(£l,U,V) 
such that 

d h (H*j,T$)<j( 7u h 2 + 7p h). (3.3) 

Proof. We have 

Tp - H h = H 2h ir 2h + (3(1 - w 2h ) - H$ 

= ((31 + lu Gf + lp G 2 p h )ir 2h + (3(1 - n 2h ) - ((31 + lu G h u + lp G h p ) (3.4) 
= lu(Gl h TT 2 h - G u ) + i P (G 2 p h TT 2h - Gp). 

We write 

G 2h n 2h - G^ = (G 2 t h — G u )n 2h + (G u — G^tt^ + G u (w2h - I)- 
For any g e Vh we have 

|((Gf - G u )g,g)\ = \(W 2h U 2h - WU)g J)\ = \\\U 2h g\\ 2 - \\Ug\\ 2 \ 

<\\(U 2h -U)g\\(\\U 2h g\\ + \\Ug\\) < Ch 2 \\g\\ 2 , 

which implies 

\\(Gl h ~G u )g\\<Ch 2 \\g\\ 
since G„ — G u is symmetric on Vh- In particular, 

\\(G 2 u h - G a )n 2h f\\ < Ch 2 \\f\\ V/e V h . 
Similarly, it can be shown that 

II(g„-g£)W1| <ch 2 \\f\\. 

Finally, we estimate 

\\G h u (I - n 2h )f\\ = \\W h U h (I ~ -K 2h )f\\ < C\\U h (I - n 2h )f\\ 

< C\\U(I - n 2h )f\\ + C\\(U-U h )(I - 7T 2h )f\\ 

J27T3J , (I2T9J1 _ _ 

< C 1 h 2 \\f\\+C 2 h 2 \\(I-7r 2h )f\\<Ch 2 \\f\\, V/ G Vh. 

Combining the last three estimates, we obtain 

\\(Gl h n 2h - G u )f\\ < Ch 2 \\f\\ V/e Vh. 
The second term in Q3.4j) is estimated in a similar way, to obtain 

|| (Gf n 2h - G p )f\\ < Ch\\f\\ V/ e V h , 
which together with the previous estimate yields 

\\(T% - H%)f\\ < C( lu h 2 + lp h)\\f\\ V/e V h . 
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It follows that 



(H h 8 fJ) 



< C 



7.(jJi/,/)+7,(ff;/,fl+/?ii/n 2 " f 



Assuming Cf3 1 (j u h.Q + J p h ) = a < 1, and < /i < h we conclude that T% is 
symmetric positive definite, and 



sup 

fev h \{o} 



In 



(H h 8 f,f) 



^ |ln(l-a)| 

< sup 

a /ev h \{0} 



{H h JJ) 



< 



\hx{l-a)\C. l2 

a p 
C 

where we also used that for a € (0, 1), x G [1 — a, 1 + a] we have 

ln(l + a), ,, llnfl - a)| , 
v - \l-x\ < [ In £tr| < - — i -|l-a;|. 



□ 



A consequence of Theorem 13.11 that legitimizes the use of Vl as a preconditioner 
for the Hessian, is stated in the following corollary. 

Corollary 3.2. There exists a constant C = C{Q.,U,V) such that 



d h (L h ,(H%r i )<j( lu h 2 + lp h). 



(3.5) 



forh<h (p,n,U,V). 

Note that, for fixed (3, since the spectral distance between the operators decreases 
with h J. 0, the quality of the preconditioner increases with h 4, 0; this is different from 
classical multigrid preconditioned for elliptic problems, where the spectral distance 
is bounded above by a constant that is independent of h. 

Next, we briefly describe how to extend the two-grid preconditioner to a multi- 
grid preconditioner that exhibits the same optimal-order quality (|3.5[) and is less 
costly to apply. If we use the classical V-cycle idea to define recursively a multigrid 
preconditioner, the resulting preconditioner is suboptimal, that is, the quality of the 
preconditioner does not improve as h 4, 0, it is simply mesh-independent. To construct 
an improved preconditioner we introduce the operators 



and 



Qh ■ £(V 2h ) -)■ C(V h ), G h (T) = Tn 2h + - 7r 2h ) 



Af h :C(V h )^£(V h ), M h {X)=2X-XH h gX, 



Note that Nh is related to the Newton iterator for the equation X 1 — Ha = 0, i.e., 
Afh(Xo) is the first Newton iterate starting at X$. Thus, if Xq is a good approximation 

than X Q . We follow [5] and 



for {H p 



h\-l 



then Afh(Xo) is significantly closer to (ifg) 



define the multigrid preconditioner Kg using the following algorithm: 
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1. if coarsest level 

2. else 

if intermediate level 

K h p =M h {g h {Kf)) 

else % finest level 

K h p =g h {Kf) 

end if 

3. end if 

Since the application of Afh requires a matrix- vector multiplication by TJg, which for 
large scale problems is expected to be very costly at the finest level, we prefer that no 
such matrix-vector multiplication is computed inside the preconditioner. This is the 
reason why we treat the cases of intermediate and finest resolutions differently In 
practice, neither Hg nor Kp is ever formed, so both are applied matrix-free (see [5] 
for details). 

The analysis of the multigrid preconditioner relics on the estimate for the two- 
grid preconditioner and properties of the spectral distance. We recall here Lemma 5.3 
from [5] that is needed for the analysis. 

Lemma 3.3. Let (ej)i>o and (<2j).j>o be positive numbers satisfying the recursive 
inequality 

ei+i < C(e, + a i+ i) 2 

and 

a i+ i <ai< f~ 1 a i+ i 
for some < / < 1. If ao < and if cq < ACa^, then 

e, < ACaj, Vi > 0. 

Theorem 3.4. Assume that Cf5~ x (^f u h\ + 7 P /io) < 2 -5 , where C is the constant 
from Theorem \3.1\ and Kg = (H^ )^ 1 . Then 

d h {K h p ,{H h p)- 1 )<2j{^ u h 2 + lp h), forh = 2- l h ,l>2. (3.6) 

Proof. Let hi = 2^ l h a , i = 0,...,l, and e t = e hi = d^Kp* , (i?^*) -1 ), a, = 
C/3~ 1 (7 tl /i i + jphi). We have 

d hi (Q{K^)AH h pT l ) < d hi (Q{Kf%g{{Hl h T 1 )) 

+ d hi {G{{HfT 1 UH$y 1 ) 

T d hi {G{Kf%Q{{Hf^)) + jfahl+^hi) (3.7) 

< d 2hi {Kf* , {HfT 1 ) + jiluhj + lp hi) 

< ej_i +ai, i = l,...,l, 
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where we have used the fact that d hi (Q hi (Ti),Q hi (T 2 j) < d 2hi T 2 ), for any two 
symmetric positive definitive operators Ti,T 2 G C(V 2 hi) Lemma 5.1]. 

It is shown in [5] that for any M,H £ £+(V/ H ), with dh t (M, H^ 1 ) < 0.4, we have 

d hi {M hi {M),H- 1 ) < 2d, H (M,H~ 1 ) 2 , (3.8) 

for /ij < ho- An inductive argument implies that < 0.2 for all i, provided it holds 
for e and C/3 _1 (7 u /i§ + 7p/io) < 0.1. Thus, combining (|3.8[) and (|3.7[) we obtain 

e 4 < 2(e,_i + a 4 ) 2 , i = 1, ...l-l. 

Note that the hypothesis of Lemma 3.1 are satisfied with / = 1/4 which implies 

C 2 

e, < 8^-(7„/i 2 +7 P /i l ) 2 , i = 1, 1. 
In particular, we have 

C 2 

e 2 /, <32 — ( 7u /i 2 +7 p /i) 2 . 

At the finest level, i = I, (|3.7[) becomes 

d h (K%, (H*)- 1 ) < e 2h + j( lu h 2 + 7p h), 

which combined with the above estimate yields the assertion of the theorem. □ 

Finally, let us note that an important difference between our method and the 
classical multigrid is that here the base case needs to be chosen sufficiently fine, 
whereas in classical multigrid it can be as coarse as possible, in general. 

4. Numerical results. In this section we present some numerical experiments 
that illustrate the application of the preconditioner introduced in Section [3l 

Let fl = (0, l) 2 and consider an optimal control problem of the form (|1.1[) . We 
consider a family of uniform rectangular grids with mesh size h and discretize the 
problem using Taylor-Hood Q2 — Qi elements for velocity-pressure and Q2 elements 
for the control. The problem was solved using Matlab R2010A. We perform three 
types of experiments: velocity control only (7^ = 1,7 P = 0), mixed velocity-pressure 
control (7 U = 1,7 P ^ 0), and pressure control only (7^ = 0,7 P = 1). 

First, we summarize the numerical results obtained for "in- vitro experiments" . As 
mentioned earlier, the Hessian matrix is dense, therefore it is never formed in practice, 
and the matrix-vector products in the preconditioned conjugate gradient (PCG) are 
implemented matrix-free. However, in order to evaluate directly the spectral distance 
between the Hessian and the proposed two-level preconditioner, we formed the matri- 
ces for moderate values of the mesh size h. In Table |4~T1 we present the joint spectrum 
analysis for (3 = 1, 7,, = 1, and j p = with dh = max{|lno;| : a € <j(Hp , Tg)}, 
where cr(A, B) denotes the set of generalized eigenvalues of A,B. The results indicate 
optimal third-order convergence. In Table l4~2l we present similar results for the case 
of pressure control only. In this case we observe an optimal quadratic convergence 
rate. We note that our computational results show a better behavior than predicted 
by Theorem 13. II We think this is due to the particular type of convex domain chosen 
here, for which the Stokes problem has better regularity than what was assumed in 
Theorem 13.11 and also to the use of quadratic elements for the control. 
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Table 4.1 

Joint spectrum analysis for /3 = 1: velocity control only (■y u = 1, 7 P = 0). 



h 


dh 


dihjdh 




1.0274 x 1(T 4 


N/A 




1.3308 x 10- 6 


7.7205 


2 _4 


1.3883 x 10- B 


9.5858 


2- b 


1.5834 x 10~ y 


8.7675 



We also remark that the numerical estimates of dh in case of pressure control 
only were obtained with a discretization that represents faithfully the finite element 
formulation in Section 12.11 However, in practice, instead of using average-zero pres- 
sures it is convenient to use a pressure space where the pressures are set to zero at 
a fixed location, e.g., a corner. If we compute the operators Hp, Tp using the latter 
space, we note that <r(Hg,Tp) contains exactly two generalized eigenvalues that are 
of size 0(1), instead of O(h), as predicted by Theorem 13. II If o{H^, Tt) is the subset 
of &{Hp, Tg ) obtained after excluding the two generalized eigenvalues that are 0(1), 
and dh = max{|lna| : a e a(Hp,Tp)}, then dh = 0(h), as the theory predicts for 
dh- 

Table 4.2 

Joint spectrum analysis for /3 = 1: pressure control only (■ju = 0, 7 P = 1). 



h 


dh 


d%hldh 




1.9613 x 10~ 2 


N/A 


2" 3 


5.6686 x 10~ :i 


3.4599 


2 _4 


1.6703 x 10- ;i 


3.3937 


2- b 


4.3735 x 10~ 4 


3.8192 



The next type of numerical experiments regard the solution of the control problem 
(jl.ip . Specifically, we compare the number of iterations required to solve the linear 
system (12.61) with unpreconditioned CG to the case when CG is used with a multilevel 
preconditioner with 1—4 levels (depending on resolution, see more comments below). 
For the results presented here, we chose the target velocity 

u d - {-2x 2 y{l - xf{l - 3y + 2y 2 ),2xy 2 {\ - y) 2 {\ - 3x + 2x 2 )) T 

and the target pressure 

Pd = cos 7tx cos iry . 

In Table 14.31 we summarize the results obtained for velocity control only (<y u = 1, 
7 P = 0), mixed velocity-pressure control (j u = 1, j p = 10~ 5 , 10~ 4 , 10~ 3 ), and pressure 
control only (fy u = 0, j p = 1); for each case a range of values for j3 is chosen. 

The choice of the values for 7 P in the mixed velocity-pressure control is justified 
by the data in Table l4~4l where we show for each case (and fixed resolution h = 1/32) 
the relative error in the recovered data for the velocity E$ — Hw 1 ,™ 111 — u<j||/||ud|| and 
pressure E p = ||p™ m — Pd||/|bd||- As can be seen from Table l4~4l for pure velocity 
control, the pressure is not recovered at all (E p ~ 1), while for pure pressure control 
the velocity is not recovered (Es ~ 1). As expected, for mixed control both pressure 
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and velocity are being recovered in the sense that both and E p decrease with 
P 4- 0. However, for 7 p = 10~ 5 the relative velocity error E$ is one order of magnitude 
smaller than E p , so velocity is better recovered than pressure. For 7 P = 10~ 4 , Eq and 
E p are of comparable size (within a factor of 2), and if 7 P = 10 -3 , then the situation 
is reversed with E p being one order of magnitude smaller than Eg. 

We now return to the actual results in Table 14.31 First note that for all cases 
unpreconditioned CG solves the system (12 .6[) in a number of iterations that appears 
to be almost mesh-independent, and is clearly bounded with h 4- 0. While this is 
certainly consistent with our remark from Section [5] regarding the mesh independence 
of the condition number of Hp, the relatively low number of CG iterations is due to the 
fact that the continuous counterpart of Hp is compact, which implies that the number 
of eigenvalues of Ha which are away from /3 is small. Also in accordance with (|2.18[) . 
the number of unpreconditioned CG iterations increases with j3 J, 0, which is expected 
since the number of relevant eigenvalues of Hp increases as (3 J, 0. For velocity control 
only (the top part of Table I4.3[) we observe a significant reduction in the number 
of iterations when multilevel preconditioncrs are used, as well as a decrease in the 
number of iterations with mesh size. For example, for f3 = 10~ 7 and h = 2 -8 , the 
number of preconditioned iterations with a four-level preconditioner is significantly 
smaller than in the unpreconditioned case, i.e, 3 iterations vs. 78 iterations. Although 
a preconditioned iteration is more expensive than an unpreconditioned one, for large 
problems the overall cost of the preconditioned solver is much lower than of the 
unpreconditioned one, as can be seen from Table 14.51 At the other end of these 
experiments (bottom of Table I4.3[) we see the cases of pressure control only. We also 
see the decrease in number of iterations with h J, when comparing the behavior of, 
say, 3-grid preconditioners at different resolutions: for example, for j u = 0, 7 P = 1, /? = 
10 -3 the number of 3-level preconditioners dropped from 14 (h = 2~ 6 ) to 10 (h — 2~ 7 ); 
similar results are consistently observed throughout Table 14.31 However, for the case 
of pressure control only, the drop in number of iterations from unpreconditioned CG 
to multilevel preconditioned CG is not as dramatic as for velocity control only. As can 
be inferred from Table l4~3l the efficiency of the multilevel preconditioned CG versus 
unpreconditioned CG measured as the ratio of the number of iterations gradually 
decreases with the increase of the ratio 7 P /7 U (for otherwise comparable experiments), 
as predicted by Theorems 13.11 and 13.41 

In our implementation we have used direct methods for solving the Stokes system, 
and we actually constructed the base-case Hessian and stored its inverse. These 
choices have limited our computations to h — 2~ 8 and /ibasc = 2~ 5 , which is why 
for h = 2~ 7 we were unable to test the two-grid preconditioner (this would have 
required saving a base-case Hessian for h = 2~ 6 ), and for h = 2~ 8 we were only able 
to compute using a four-level preconditioner. While we already commented on the 
positive side of the numerical results, it is worth noting the pitfalls: if the coarsest 
grid is too coarse, then the quality of the preconditioner declines to the point that it 
is hurting the computation. This can be seen in the groups of columns and rows in 
Table l4~3l corresponding to h — 2~ 6 and /3 — 1CP 6 , 10~ 7 : the use of too many levels 
results in a spike in the number of iterations to the point of non-convergence (within 
a maximum number of 100 iterations allowed). 

Finally, we would like to comment on the robustness of our algorithm with respect 
to the accuracy of the Stokes solve. For large-scale problems the Stokes system on 
the finer grids is expected to be solved using iterative rather than direct methods, 
which reduces the accuracy of computing matrix- vector multiplications for the Hessian 
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matrix Hg. We have repeated our numerical experiments in that, except for the 
coarsest scale, we replaced the direct solve of the Stokes systems with a preconditioned 
MINRES solve, as described in [7J. For the case of pure velocity control (jp = 0) we 
found no significant change in the number of iterations in Table 14.31 However, for 
the case of mixed- and pure pressure control (j p ^ 0) the quality of our algorithm 
appeared to decline significantly. We identified as the primary cause for this behavior 
the fact that, even when the velocity variables are well resolved, that is, the relative 
error between the solutions obtained via direct vs. iterative methods is on the order 
of 10 -8 , the relative error in the pressure terms can be quite high (10 -2 -10~ 4 ). Since 
our algorithm relies on the ability to compute the operator Vh with sufficient accuracy, 
we are not able at this point to draw conclusions with respect to the influence of using 
iterative methods on our algorithm for mixed- and pure pressure control. 

Table 4.3 

Iteration count for multilevel preconditioners; "nc" means "not- converged" . The tolerance is 
set at Id" 12 . 



h 




2- 


6 




2" 7 


2~ 


s 


num. levels 


l 


2 


3 


4 


1 


3 


4 


1 


4 






7 


u = 1, 7 P 


= 










P = lfj- 4 


7 


3 


3 


3 


7 


2 


2 


7 


2 


13 = 1Q- 5 


13 


3 


3 


3 


13 


3 


3 


14 


2 


13 = 1(T 6 


29 


4 


4 


6 


29 


3 


3 


32 


3 


p = 1(T 7 


74 


5 


7 


62 


75 


3 


5 


78 


3 






7« 


= 1, 


1p = 


10~ 5 










13 = icr 4 


10 


4 


5 


5 


10 


5 


5 


11 


5 


P = l(T 5 


20 


5 


6 


8 


20 


6 


8 


22 


7 


p = icr 6 


45 


7 


8 


13 


44 


7 


9 


45 


9 


p = lO" 7 


112 


9 


14 


nc 


113 


9 


14 


119 


11 






lu 


= 1, 


1p = 


io- 4 










p = 1(T 4 


10 


5 


6 


8 


11 


5 


7 


11 


7 


p = 10~ 5 


21 


7 


7 


9 


23 


7 


9 


24 


9 


p = io- 6 


48 


8 


10 


16 


48 


9 


13 


49 


11 


p = io- ? 


122 


13 


17 


nc 


126 


11 


22 


133 


20 






lu 


= 1, 


1p = 


10- 3 










p = icr 4 


12 


6 


7 


9 


12 


7 


9 


13 


9 


P = l(T 5 


25 


8 


9 


13 


26 


9 


13 


28 


11 


p = io- 6 


61 


11 


17 


33 


63 


12 


22 


65 


20 






7 


u = 0, 7 P 


= 1 










p = 10" 1 


8 


5 


7 


9 


8 


6 


8 


8 


8 


p = io- 2 


13 


7 


9 


12 


13 


8 


11 


13 


11 


p = icr 3 


27 


10 


14 


nc 


26 


10 


15 


27 


14 



5. Conclusions. In this article we introduced Schur-based two- and multigrid 
preconditioners for the KKT system associated to a Tikhonov-regularized optimal con- 
trol problem constrained by the Stokes system. We showed that, if the Stokes-system 
is discretized using a stable pair of finite elements, the preconditioner approximates 
the reduced Hessian of KKT system to optimal order with respect to the convergence 
order of the finite element method. As a consequence, the number of preconditioned 
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Table 4.4 

Relative error for recovered data for velocity control Eg = |m™ u1 — "d||/||"d|| and pressure 
E P = \K ln -Pd\\/\\Pd\\- 



In 


7 P 


Eu 


E p 


Eg 


E p 


En 






10 


-5 


10 


-a 


10 


-7 


1 





4.07e-2 


w 1 


8.23e-3 


w 1 


1.88e-3 


« 1 


1 


10~ 5 


4.32e-2 


4.32e-l 


1.65e-2 


2.80e-l 


8.28e-3 


7.05e-2 


1 


1(T 4 


6.66e-2 


2.75e-l 


3.30e-2 


6.85e-2 


9.55e-3 


8.44e-3 


1 


KT 3 


1.23e-l 


6.56e-2 


3.78e-2 


8.18e-3 


9.70e-3 


8.98e-4 




10 


-i 


10 




10 


-3 





1 


1.02 


2.66e-l 


1.07 


6.12c-2 


1.09 


7.34e-3 



Table 4.5 

Velocity control only: time comparison for h = 2~ 8 . No. state variables: 588290, no. control 
variables: 5222^2. 



no. levels 


1 


4 


/3 = 10~ b 

p = io- 7 


3613 s 
8953 s 


493 s + 1486 s for base case 
575 s + 1492 s for base case 



CG iterations needed for solving the optimal control problem to a given tolerance 
decreases with increasing resolution, asymptoting to just one iteration as h J, 0. The 
problem discussed in this article forms an important stepping stone towards finding 
highly efficient methods for solving large-scale optimal control problems constrained 
by the Navier-Stokes system, which are the focus of our current research. 

Appendix A. Some facts about spectral distance. 

Let (X, (•, •}) be a real finite dimensional Hilbert space and denote the complex- 
ification of X by X c = {u + iv : u,v G X}. Let C+(X) = {T G C(X) : (Tu,u) > 
0, Vu6l \ {0}} and define the spectral distance between S,T G C + {X) to be 



d x (S,T)= sup 

«i6XC\{0} 



In 



(S c w, w) 



(T c w, w) 



where T c (u+iv) = T(u)+\T(v) is the complexification of T. The following inequalities 
were proved in [5j Lemma 3.2]: 

Lemma A.l. If a G (0, 1) and z € B a (l), then 



lnfl + a), |ln(l-a)l. 
v ' \1 -z\<\laz\< - — -\1 - z\ 



(A.l) 



For \lnz\ < 5 we have 



lnzl < II - z\ < 



In; 



(A.2) 



Lemma A. 2. Let L,H G C + {X) such that 



mm{d x {L~ 1 ,H),d x {L,H' 1 )) <S 
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Then 

e s - 1 

p{I-LH) < — — min (d x (L-\H),d x (L,H- 1 )) . (A.3) 

Proof. If A 6 a (I — LH) then there exists a unit vector u 6 X c such that 
(/ — LH)u = Xu, therefore 

(1 - X)u = LHu . (A.4) 
After left- multiplying with L^ 1 and taking the inner product with u we obtain 

(Hu, u) 

(1 — A) (L~ 1 u,u) = (Hu,u) , therefore A = 1 — ' — r . 

(L L u,u) 

If we substitute v = H~ 1 u in (|A.4[) and take the inner product with v we have 

r-l„._r„. ^„^„„„x_i ( Lv > v ) 



(1 - X)H~ L v = Lv , therefore A = 1 - 



Hence, if dxiL^ 1 , H) < 5, then 

p(I-LH) < sup{|l-2| : z = (Hu, u) / {^L^ 1 u, for some u G X C \ {0}} 

£0 e 5 — 1 
< ^—jr—dx(L ,H) . 

Instead, if dx(L, H^ 1 ) < 5, then 



p(I — LH) < sup{|l — z\ : z — (Lu,u } / {H^u, u) for some u E X C \ {0}} 

< 

i 

which proves (IA.3I) . □ 



El e 5 - 1 

< e - 1 -d x {L,H- 1 ) 
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